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We study the corrections to scaling for the mass of the watershed, the bridge line, and the optimal 
path crack in two and three dimensions. We disclose that these models have numerically equivalent 
fractal dimensions and leading correction-to-scaling exponents. We conjecture all three models to 
possess the same fractal dimension, namely, d/ = 1.2168 ± 0.0005 in 2D and d/ = 2.487 ± 0.003 in 
3D, and the same exponent of the leading correction, Q, = 0.9 ± 0.1 and Q, — 1.0 ± 0.1, respectively. 
The close relations between watersheds, optimal path cracks in the strong disorder limit, and bridge 
lines are further supported by either heuristic or exact arguments. 
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I. INTRODUCTION 

The watershed, defined as the line separating adjacent 
drainage basins (catchments) , plays a fundamental role in 
water management [TH3] , landslides gHZI ; ^^'^ flood pre- 
vention [7l-(9]. From observations of watersheds in nature, 
claims about their fractality have been made already long 
ago [TU]. More recently, watersheds were investigated in 
Refs. |llH13j where their self-similarity was shown nu- 
merically for both natural and artificial landscapes. 

A fractal dimension consistent with the one of water- 
sheds was also found for optimal path cracks in the limit 
of strong disorder. Optimal path cracking has been in- 
troduced by Andrade et al. [TUfTS] as a model for the 
evolution of successive optimal paths under constant fail- 
ure. It describes, e.g., the breakdown of electrical or fluid 
flow through random media and has important applica- 
tions also in other fields of science and technology, such as 
human transportation, fracture mechanics, or polymers 
in random environments, where finding the optimal path 
is a challenge [T7H27] . 

The last of the three problems mentioned in the ti- 
tle, related with ranked percolation (RP) was recently 
introduced by Schrenk et al. |28] as a model where the 
creation of a spanning cluster is systematically delayed. 
They found that the set of "bridge bonds" (i.e. bonds 
that finally lead to spanning clusters) has a fractal dimen- 
sion very close to that of watersheds and of the optimal 
path cracks in strong disorder (see Fig. [T]). 

The appearance of the same fractal dimension in three 
seemingly very different models calls on the one hand for 
a theoretical explanation, and on the other hand for more 
precise numerical estimates. On the theoretical side, we 
might point out that the watershed (WS), the optimal 
path crack in strong disorder (OPC), and the bridge line 
(BL) in RP are all sets of sites or bonds that split the sys- 
tem into two distinct parts and seem conceptually related 
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(although not identical) to classical percolation. Yet, de- 
spite these similarities and the broad relevance of the 
models, no detailed studies of the relation between them 
are available. 

Finally we should mention that also relations to other 
physical models have been proposed, such as optimal 
paths jTBl l?5H3^ , the shortest path in loopless invasion 
percolation |30j . the infinite cluster in multiple invasion 
percolation [33] , and the surface of the infinite cluster in 
explosive percolation jMll35j. 

In this paper we explore the relation between the main 
crack (MC) of the optimal path crack in strong disor- 
der pAro] and the bridge line of RP ^. But we shall 
also explore the relations between several definitions of 
watersheds pTHT5] , since the exact definition of a water- 
shed is not unique, and different definitions turn out to 
be closely related to different subsets of the other three 
problems. We present improved estimates of the frac- 




FIG. 1. (color online) Mass M of the watershed (WS 
site/bond), the main crack (MC), and the bridge line (BL) 
as a function of the system size A*', defined as the number of 
sites (bonds) in the system, in both two and three dimensions. 
The error bars are much smaller than the symbols. The lines 
show the fractal dimensions obtained in this work. 



TABLE I. Number of samples used to obtain the average mass of the bridge hue (BL), the watershed (WS sites/bonds), and the 
main crack (MC) for different system sizes L in two- and three-dimensional systems. For the numerical analysis of corrections 
to scaling it is important to use high-precision data. Therefore, we focused on obtaining best possible statistics for the lattice 
sizes listed here, instead of increasing the number of different system sizes. 
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tal dimensions, made possible by studying in detail the 
corrections to scaling for two- and three-dimensional sys- 
tems with uncorrelated disorder. Due to the numerical 
difficulty in obtaining sufficient statistics, we omit a dis- 
cussion of the surface of the infinite cluster in discontin- 
uous (explosive) percolation [MUSS]. For all models, the 
fractal dimension df is defined through the scaling of the 
mass M, corresponding to the number of sites or bonds 
in the object, with the linear system size L, 



M - L"' 



(1) 



Due to the finite system size, corrections to scaling arise 
[36H38J that may mask the true asymptotic behavior. 
Hence, the estimated df can be improved by describing 
the size dependence of the mass as 



Mj=L'^fCT 



(2) 



where the general form for the corrections to scaling Cl 
is 

Cl = aoo + aoi^"^ + 002^^^ + ao3-^"^ + • ■ • 

+aiiL-"' + ai2L-"'-i + 013^"^^'"' + ■ • ■ 
+a2iL-^' + 022^""^"' + a23L-^^'-^ + ... 
+a„iL-^" + . . . , (3) 

with non- universal coefficients (a.y). The exponents ful- 
fill Oi < CI2 < ... < r2„ and are non-analytic (non- 
integer). They are usually independent on the geome- 
try of the lattice and only depend on the dimensionality 
|36|, I37j . Finding the same non-analytic corrections-to- 
scaling exponents for all three models will give another 
hint for the close relation between them. But, in gen- 
eral, the precise estimation of corrections to scaling is a 
difficult task. Numerical studies typically measure the 
leading correction exponent, a sub-leading correction ex- 
ponent, or an effective exponent arising from the sum of 



two or more correction-to-scaling terms '35|. Hence, a 
reliable estimate of the leading correction exponent de- 
pends on both the method and the precision of the data. 
Since in practice it is not reasonable to attempt a fitting 
with many terms of the form shown in Eq. (13 ) , we trun- 
cate the sum of correction terms as discussed in detail 
below. We first have a look at the functional form of 
the corrections to scaling that can be considered for the 
individual models given the available statistics, using a 
simple fitting and checking which amplitudes in Eq. ([3]) 
are small. Using this and truncating terms with an ex- 
ponent > 3, we define our effective corrections-to-scaling 
ansatz. By defining a fit quality, we identify the leading 
correction exponent (highest maximum of the quality) 
and obtain a highly accurate estimate for the fractal di- 
mension df . The largely improved estimate of df is the 
main focus of our numerical study, rather than obtain- 
ing the corrections with precision. We cross check the 
obtained results with a careful analysis of the local loga- 
rithmic slopes as suggested by Ziff [301 SI] ■ This method 
uses the fact that for large enough system sizes the higher 
order terms are negligible, such that the local logarithmic 
slope of the corrections to scaling should converge to the 
leading correction exponent. 

The paper is organized as follows. In Sec. [H] we de- 
scribe the models. Section UlTl introduces the corrections 
to scaling and summarizes the obtained results. The re- 
lations between the models are discussed in Sec. IIVI and 
conclusions are drawn in Sec. IVl 



II. MODELS 

In the following, we give a brief overview of the wa- 
tershed (WS), optimal path cracking (OPC), and ranked 
percolation (RP), focusing on the role of percolation in 



the numeric procedures used to determine the watershed, 
the main crack (MC), and the bridge hne (BL). For sim- 
phcity, the description is given for two-dimensional sys- 
tems (square lattices), where they lead to lines. The 
extension of the discussed models to higher dimensions, 
where they lead to (hyper) surfaces, is straightforward 
and has been done in Refs. [131 [T^[25] . 



A. Watersheds 

Watersheds are the lines separating adjacent drainage 
basins and play a fundamental role in many fields [1|- 
[Sj. Although the intuitive notion of a watershed seems 
obvious, the choice of a precise definition is rather subtle. 
Indeed, in the previous literature (see |42| for a review) 
several definitions have been used, none of which seems 
optimal. Moreover, as we shall see, the choice of the most 
efficient algorithm for simulating a watershed depends on 
the precise definition, and different definitions - although 
corresponding to the same "macroscopic" objects - are 
more or less directly related to the other two problems 
discussed in this paper. 

Following TTHT5] , we shall discuss in the present paper 
two main definitions, the bond model and the site model, 
and in addition a variant of the latter, the great wall 
m,odel (called flooding method in [11'). As we shall see, 
the natural algorithm for the bond model is one where we 
follow the run-off from top to bottom, while the natural 
algorithms for the site models 'flood' the catchment areas 
from their outlets to the top. 

We consider uncorrelated artificial landscapes mapped 
on a square lattice of size L x L as digital elevation maps, 
where each site i = {x,y) represents a small square area. 
The height hi at each site i is drawn randomly from a 
common distribution in such a way that hi > 0. The pre- 
cise form of the distribution is irrelevant, provided it is 
continuous so that, with probability one, hj ^ hi, ^j^i- 
Boundary conditions are periodic in the horizontal direc- 
tion, but free vertically. Thus water can run across the 
lateral sides in both directions (depending on which of 
the neighboring sites is higher), while it can only flow 
outwards from the top {y = L — 1) and bottom {y = 0). 
The latter could be modeled more explicitly by adding 
two more rows (with y = L and y ~ —1) where all sites 
have height h = 0, and which act as sinks. The parts of 
the landscape that drain to either of these two sinks are 
their catchment basins, while the line separating the two 
catchment basins is the watershed. 

Water flows always from a higher site to a lower one, 
but the bond and site models correspond to different as- 
sumptions how this happens in detail. In the bond model, 
water flows from any site only to its lowest neighbor, while 
it flows to all neighbors in the site model. Thus, each site 
belongs in the bond model to a unique catchment area, 
and the watershed must be formed by bonds of the dual 
lattice which cut bonds that join sites in different catch- 
ment basins. It is easy to see that a watershed defined 



this way must be a single connected and unbranched path 
that has no loops except for the fact it is periodic in the 
horizontal direction (and is thus one big loop). Moreover, 
determining the catchment basin for any site is trivial: 
one just has to follow the unique run-off path. 

In contrast, sites do not have unique run-off paths in 
the site model. Let us call a site with more than one 
lower neighbor a diversion site. At each diversion site, 
the run-off path branches, so that the total run-off pat- 
tern of any site is a tree. Moreover, branches of this 
tree might end in both sinks, in which case the site can- 
not be in either catchment basin. Such sites must be- 
long thus to the watershed, while sites which drain into 
one unique sink form the catchment basins. Finally, two 
adjacent sites i and j cannot be in different catchment 
basins (since either hi < hj or hi > hj). Therefore the 
entire watershed must be formed by a single loopless and 
unbranched chain of sites, that is connected in the sense 
that adjacent sites must be either nearest or next-nearest 
neighbors. 

While it is in principle possible to follow the entire run- 
off trees in case of the site model, it is not very practical 
and easy. Thus it is more efficient to determine the water- 
shed by a flooding algorithm, where the catchment areas 
are determined by moving inward & upward from the 
sinks. Below we shall describe two such algorithms that 
differ in details. On the other hand, for the bond model 
it is very efficient and easy to follow the run-off, as de- 
scribed in [TT]. We first determine the catchment basins 
for the sites on a search line {x = 0, y) with y = 0, 1, 2 . . .. 
The first ones will drain to the bottom. After we have 
found the first site draining to the top, we have also the 
first bond in the watershed. Starting from this bond we 
can then construct the entire watershed recursively, by 
following the run-off paths from the sites adjacent to one 
of its endpoints. 

For the site model we flood simultaneously two inva- 
sion percolation clusters growing inward from the top and 
bottom rows. Let us call Bt{h) and Bi,{h) the boundaries 
of these clusters, when the flood has height h. More pre- 
cisely, Bt{h) {Bb{h)) is the set of all sites i with height 
hi > h, and with at least one neighbor j having hj < h 
and being in the top (bottom) cluster. Starting with 
h = 0, we increase h continuously, each time incorpo- 
rating a boundary site into the corresponding cluster, as 
soon as it gets flooded - provided this site does not belong 
to both boundaries. A site belonging to both boundaries 
obviously drains into both basins and is thus part of the 
watershed. 

When reaching the first site on the watershed, we have 
two options. In one, we flood it like any other site, but 
take care that any site draining into it must also be in the 
watershed. Thus, when increasing h further, we have to 
distinguish between sites that get flooded from neighbors 
that all belong to the top basin, sites that get flooded 
only from neighbors that all belong to the bottom basin, 
and sites that get flooded either from both or from a site 
in the watershed. The flrst belong to the top basin, the 



second to the bottom basin, and the third to the water- 
shed. The algorithm stops when the entire landscape is 
flooded. This gives the site model proper, and is meant 
whenever we speak of the 'site model' in the following 
sections. 

Alternatively, we can prevent sites on the watershed 
from being flooded by increasing their height to a value 
larger than any other hi in the entire landscape. In this 
way the two floods are kept separated, and we can con- 
tinue flooding without any further modification. The al- 
gorithm stops when the entire landscape is flooded ex- 
cept for the watershed sites. These sites form then a 
connected wall (or dam), whence the name great wall 
model. We will not present data obtained with this algo- 
rithm directly, but it is most closely related to the models 
discussed in the next two subsections. 

The mass M of the watershed (WS) is deflned as the 
number of bonds (sites) forming the watershed. No- 
tice that we do not consider the watershed as a three- 
dimensional object (with height as third dimension), but 
as 2-dimensional, see Eq. ^. 



B. Optimal Path Crack 

The optimal path crack (OPC) was introduced by An- 
drade et al. P^l - fTB] and is obtained in the following 
way. We start with a square lattice of size L using free 
boundary conditions in the vertical direction and periodic 
boundary conditions in the horizontal one. A random en- 
ergy is assigned to each site and the energy of any path in 
the system is deflned as the sum of the energy of its sites. 
In particular, the optimal path is the one among all paths 
connecting the top and bottom boundary of the system 
with the lowest total energy. Once the first optimal path 
is determined, the site in the optimal path having the 
highest energy is identified and removed. This is equiva- 
lent to imposing an infinite energy to this site. Next, the 
optimal path is calculated among the remaining acces- 
sible sites of the lattice, from which the highest energy 
site is again removed. The process continues iteratively 
until the system is disrupted and no further path can be 
found. The set of removed sites then defines the optimal 
path crack (OPC). The OPC is dependent on the type of 
disorder, but in the limit of strong disorder, it is localized 
in a single line, denoted as the main crack (MC), with 
mass M given by the number of cracked sites. From this 
point on, we consider the OPC only in the limit of strong 
disorder and, for simplicity, just refer to it as main crack 
(MC). 

In the strong disorder limit, the model is equivalent 
to the great wall model, with h corresponding to the 
random energy and the main crack corresponding to the 
great wall. 



C. Ranked Percolation 

Ranked percolation is a new percolation model intro- 
duced by Schrenk et al. 28] in which the creation of a 
spanning cluster is suppressed. In this model bonds or 
sites are occupied randomly, except for bridges that are 
bonds/sites which, when occupied, would create a span- 
ning cluster, i.e. a cluster connecting top and bottom 
edges of the system. In the following, we focus solely on 
the case where bridges are never occupied (in the more 
general model of [28] they have a probability pi, of being 
occupied that is smaller than the probability for other 
bonds/sites; in this notation, the present simulations cor- 
respond to pb = 0). While the original studies were done 
for bond percolation, we consider here site percolation. 
Similarly as in the bond case, we start with an empty 
square lattice of size L x L, choose sites uniformly at 
random and occupy them. If two neighboring sites are 
occupied, they are considered to be connected and to 
belong to the same cluster. In contrast to standard site 
percolation, whenever the occupation of a site would lead 
to a spanning cluster, this bridge site is blocked. The pro- 
cess proceeds until all sites are occupied or blocked and 
the system is disrupted into two parts. The separating 
bridge line (BL) is formed by the set of bridge sites. 

Cieplak, Maritan, and Banavar [29] have studied this 
line in a different context and argued that the occupa- 
tion procedure is equivalent to the following: Randomly 
assign an energy to each site, rank order them by increas- 
ing energy, and occupy them according to their rank - 
except when the site to be occupied is a bridge site. In 
that case the site is not occupied ever. Seen this way, it 
transpires that also ranked percolation is equivalent to 
the great wall model, except for the fact that sites are 
'fiooded' in different order and the algorithms suggested 
by the two models are very different. 

Finally, let us point out that the bond version of ranked 
percolation is not strictly equivalent to the bond model 
defined in subsection |II A[ but corresponds to a bond 
model on a slightly different lattice [28] . 



III. CORRECTIONS TO SCALING 

We perform extensive numerical simulations of the de- 
scribed models measuring the mass M of the watershed 
(WS), the bridge line (BL), and the main crack (MC) for 
different (linear) system sizes L. For details about the 
considered system sizes and the corresponding number 
of samples, see Tab.|lj The obtained masses are shown in 
Fig. [l] as a function of the system size N, namely N = L"^ 
for sites and N = dL'^ - {2d - 1)L'^-^ for bonds (the 
second term arises due to the solid walls in the vertical 
direction), where d is the dimensionality of the system. 
Although this is not visible in Fig.[l] the masses of the 
BL and MC are equal within the error bars. Those of WS 
site and WS bond are different from the masses of BL and 
MC. Nevertheless, we observe all of them to follow very 
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FIG. 2. (color online) Corrections to scaling Cl = M/N'^f^'^ 
of the watershed (WS site/bond), the main crack (MC), and 
the bridge line (BL) as a function of the system size A'^, defined 
as the number of sites (bonds) in the system, in two dimen- 
sions. The fractal dimension df ~ 1.217, consistent with the 
more precise estimate obtained later, has been chosen such 
that Cl converges to a constant value for large N. The error 
bars are typically much smaller than the symbols. The lines 
show fits of truncated versions of Eq. pl to the data, which 
is divided here by aoo to show the matching of the scaling 
behavior of the different models for large A*'. 



similar scaling behaviors. The true asymptotic behavior 
for the mass scaling, Eq. ^ , is masked by corrections to 
scahng arising due to finite system size [36l[37]. Hence, 
the estimate of the fractal dimension dj can be improved 
by considering these corrections explicitly, see Eq. ([3]) . In 
the following, we first analyze the general ansatz to find 
the number of distinguishable correction exponents and 
if there are vanishingly small amplitudes. This results 
in simplified functional descriptions of the corrections to 
scaling in 2D and 3D, which are then studied by two 
different techniques in order to obtain highly accurate 
estimates of the exponents. 



A. Ansatz for Corrections to Scaling 

To understand the structure of the data, we study 
least-square fits of different truncated versions of Eq. ([S]) 
to the corrections to scaling Cl — MjU^s in 2D, where 
df ~ 1.217 has been chosen such that Cl converges to 
a constant value for large L (see Fig. [2]) . This choice of 
dj is consistent with the more precise estimates obtained 
later. Using different numbers of exponents r2„ and vary- 
ing numbers of expansions, we find for all models that, 
with the current precision, we cannot resolve correction 
terms of an order higher than XjL? . In the following, we 
therefore truncate the expansions by setting aij = 0Vj>2. 
For the case of WS site, we obtain reasonable fits down 
to fairly small L using a set of two exponents (n = 2), 
yielding Jli « 0.6 and 1^2 ~ 0.9, while ai2 seems to be 
small and also the amplitudes of the analytic terms seem 



to be small and unresolvable (aoi ~ 0, ao2 ~ 0). It 
is important to note that, despite these findings, U,2 is 
still compatible with unity. For WS bond, MC, and BL 
we obtain similar results, although an seems to be very 
small in all the three models. As shown in Fig. [2] our 
fits match Cl fairly well for the models. Hence, defining 
a; = r^i and fJ = ri2 the (visible) corrections reduce to 
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with Lu « 0.6 and 51 w 0.9, while an is large only for 
the WS site model. The latter fact will be discussed 
in section ITVl We note that we did not find evidence of 
logarithmic corrections. Figure [2] showing the rescaled 
data, confirms that the the corrections considered here 
capture the behavior of the data. 

In 3D, we find by a similar study, that the corrections 
to scaling of all four models can reasonably well be de- 
scribed by a single correction term such that we can write 



C'l — 0,00 + aiiL 



(5) 



with Q ss 0.9, but compatible with unity. A simple least- 
squares fit of the ansatz given by Eq. (|4| (Eq. ^ in 3D), 
to the data to obtain the coefficients, df, fi, and/or lu 
directly can be ambiguous. Dependent on the choice of 
the initial values for the fit parameters (the coefficients 
and exponents), a fit could even lead to an estimate of 
51 or a; reflecting higher order corrections instead of the 
leading ones. To overcome this and improve the accu- 
racy, we discuss, in the following, a method that explores 
the parameter space by varying the exponents in a given 
range and analyzing the quality of the corresponding fits. 
If one would attempt to fit an ansatz containing at the 
same time terms with variable exponents and analytic 
corrections to the data [formally similar to Eq. ([3])], in- 
terference among the terms would be possible when the 
variable exponent is close to unity. The fact that 51 is 
close to unity, does not affect our procedure, since the 
corrections given in Eqs. (H| and ([5]) do not contain an- 
alytic terms explicitly. The results from this method are 
then cross checked with a second method, which allows to 
estimate the leading correction from the convergence of 
the local logarithmic slopes in the reduced mass M L~'^f . 



B. Fit Quality Method 



The output of a fit of the ansatz in Eq. (B or in Eq. (Is]) 
to the data of the reduced mass M L^'^f can be sensitive 
to the initial conditions. We, therefore, perform a more 
systematic study as follows. To have a good control over 
the actual fitting, we use Eq. (l2| , (|4]) , and (Is]) in the fol- 
lowing form. 
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FIG. 3. (color online) (a) Inverse of the quality Q as a func- 
tion of Q' for different values of a, as obtained from fits of 
the ansatz, Eq. m, to the reduced mass M L~°' for the wa- 
tershed on bonds in 2D with sizes as indicated in Tab. |T] The 
vertical lines give the position of the global minimum in 1/Q 
(solid) and the estimated error (dashed), (b) For the same 
system as in (a), the minimum value 1/Qioc as a function of 
a is shown, where Qioc is obtained from curves 1/Q{Q') for 
a given a, like those shown in (a). The vertical lines high- 
light the value of a at the global minimum 1/Qmax (solid) 
and the estimate for the error (dashed). The error bars are 
determined from the width of the minima. The vertical lines 
show the estimate df = 1.2168 ± 0.0005 for the fractal di- 
mension of the watershed on bonds and the horizontal ones 
the corresponding leading correction Q = 0.95 ± 0.05. These 
exponents were obtained from the analysis of a single model 
(WS bond). By combining the results for different models, 
we obtain more accurate estimates for the exponents. 



in 2D and 3D respectively, with fixed values of a, a;', and 
il' and estimate the Quality Q = n/x^, where n is the 
number of degrees of freedom of the fit, i.e. the number 
of system sizes used in the data (see Tab. IT]) minus the 
number of fit parameters (here the number of resolvable 
amplitudes), and x^ is the (weighted) mean square devi- 
ation of the fit. The quality Q is a function of a, w', and 
57', but, as the terms of oj only have visible amplitudes for 
the WS site model in 2D, we drop hereafter the depen- 
dence of Q on w' and fix w' = 0.6. Since uj' is fixed, only 
one single correction exponent, f2', is adjusted, avoid- 
ing fitting simultaneously multiple exponents. Now, Q 
should be maximal for a = df and n' = n, as the leading 
correction gives the dominant contribution compared to 
higher order ones. As a matter of convenience, we use the 
inverse of the quality 1/Q, which is minimal for a = df 
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FIG. 4. (color online) Inverse of the quality at the minimum 
1/Qioc as a function of a for the different models in 2D. The 
inset shows for each model the inverse of the quality Q as a 
function of Q' with a fixed to its value at the global minimum 
l/Qmax- The vertical lines show the averages df = 1.2168 ± 
0.0005 and Q = 0.9 ±0.1 of the estimates for the fractal 
dimension and for the leading correction, respectively. 



and n' = n. The procedure to obtain df and f2 for a 
given model is to measure the inverse quality 1/Q(a, $7') 
of a fit of the proper ansatz to the data. We first choose 
a value of a and then derive 1/Q as a function of the 
exponent of the leading correction, scanning in the range 
< 17' < 2 with a step size Sfl' = 0.015. The obtained 
curve, see for example Fig. Isla) for WS bond, typically 
has a (local) minimum 1/Qioc(q^): which marks the best 
fit of the leading correction r2ioc(Q;) for the chosen a. 
The error Ariioc(a) is estimated from the width of the 
minimum. In two dimensions, analyzing these minima 
l/Qioc(a) by varying a in the range 1.212 < a < 1.220 
with steps of size 5a = 0.00025, yields an estimate of the 
global minimum l/Qmax and the fractal dimension df. 
The error bar in df is also determined from the width of 
the minimum [see, e.g., Fig.lslb)]. We repeated this anal- 
ysis for the watershed on sites, the main crack, and the 
bridge line (see Fig. [4]) and the corresponding estimates 
are summarized in Tab.lTlj The obtained values all agree 
with each other within the error bars. The ones for the 
MC, due to the low statistics, seem to differ more. Nev- 
ertheless, based on the similarity of the numerical values, 
we estimate by combining the results for all models that 
in 2D df = 1.2168 ± 0.0005 and O = 0.9 ± 0.1 for all 
models. The values and error bars have to be obtained 
by a reproducible procedure. We used the intersection of 
the estimated intervals for all models (Tab. In]). The value 
obtained for fi is close to unity, which suggests that the 
leading correction (the second leading correction for WS 
site) is likely to be the analytic correction ft = 1. 

We applied a similar analysis to the data obtained 
in three-dimensional systems, scanning 51' in the range 
< 57' < 2 with a step size Sfl' = 0.015 and a in the 



TABLE II. The fractal dimension df and the exponent of the 
leading correction fl of the bridge line (BL) and the watershed 
(WS sites/bonds) for 2D and 3D, as obtained from a similar 
analysis as done in Fig. |3] for the WS bond case. The main 
crack (MC) result is only shown for 2D. 



model 


d 


df 


Q. 


WS bond 


2 


1.2168±0.0005 


0.95±0.05 


WS site 


2 


1.21705±0.00075 


0.91±0.19 


BL 


2 


1.2166±0.0015 


0.87±0.08 


MC 


2 


1.2166±0.0045 


0.86±0.11 


WS bond 


3 


2.4865±0.0025 


0.96±0.10 


WS site 


3 


2.4865±0.0025 


0.98±0.09 


BL 


3 


2.4878±0.0025 


1.06±0.16 



range 2.450 < a < 2.535 with steps of size da — 0.0025. 
As before, the detailed analysis is done like is shown in 
Fig. Isl). For the case of the main crack in 3D, no con- 
clusive results could be obtained with our method, but 
the obtained masses are within their error bars equiv- 
alent to those measured for the bridge line. We show 
in Fig. [5] only the results obtained for the watershed on 
bonds, on sites, and the bridge line. Like in 2D, the 
obtained estimates for df and fl agree within the error 
bars. Therefore, we estimate df = 2.487 ± 0.003 and 
il = 1.0 ± 0.1 for three dimensions. As in 2D, the value 
of the leading correction is likely to be the analytic one 
n = 1. Given this, for 2D and 3D, we also analyzed the 
data fixing J7 = 1. The obtained values for the fractal 
dimensions and their error bars are consistent with the 
ones reported in Tab. |lTj therefore, the possibility of fi 
being analytical cannot be discarded. 

The estimates of the fractal dimension for the different 
models are in agreement with the ones found in previous 
works for the watershed (1.211±0.001 [TT and 2.48±0.02 
[T3j), the main crack (1.215 ± 0.005 and 2.46 ± 0.05 iH- 
Hg), the bridge line (1.215±0.003 and 2.50±0.02 [MSO]), 
and the perimeter of the infinite cluster in discontinuous 
percolation (1.23±0.03 [34 and 2.5±0.2 [SSj). The value 
1.211±0.001 given in Ref. [TT] for the fractal dimension of 
the watershed in two dimension seems to underestimate 
the error bar. 



C. Local Logarithmic Slope 

Another approach to estimate the leading correction- 
to-scaling exponent fl is to calculate the local logarithmic 
slope of the reduced mass Cl{c() = MlL~°', i.e.. 



n^st[L,a) = -log2 T- 7- 

\CL/2{a) -Ci/4(a 



(7) 



Taking L relatively large such that higher order correc- 
tions are negligible, fiest converges to fi for a = df (see, 
e.g., Refs. [101111]). Due to the uncertainty AM^ in the 
average of the mass M^ , there are in the estimate of the 
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FIG. 5. (color online) Inverse of the quality at the minimum 
1/Qioc as a function of a for the watershed on bonds, the 
watershed on sites, and the bridge line in 3D. The inset shows 
for each model the inverse of the quality Q as a function of 
Q' with a fixed to its value at the global minimum l/Qmax- 
The vertical lines show the averages df = 2.487 ± 0.003 and 
57 = 1.0 ± 0.1 for the estimates for the fractal dimension and 
for the leading correction. 



local slope systematic errors of the form 



/ df^e 



-ACr 



VdCi/fe-^^/'^ 



A^UL) = J2 

fe={l,2,4} 

^ f {ACL? + {ACL/2r 

\ {Cl - Ci/2)2 

nACL/2r + {acl/4) 



V (Cl/2 - Cl/aY 



{ACl/2)' 



{Cl - Cl/2){Cl/2 ~ Ci/4) 



(8) 



where ACl — L^^AMl. We omitted here the a depen- 
dence. This error heavily depends on the precision of 
the single mass measurements and, therefore, statistics 
considerably higher than for the fit quality method are 
needed, especially for the larger system sizes. We focused 
mainly on improving the statistics for the watershed on 
bonds and for the bridge line, where larger systems can 
be addressed. 

In Figs. [6] and [7] we show ftcst with a = 1.2168 and 
2.487 for two- and three-dimensional systems, respec- 
tively. In both figures, only values of J7est with Afiost < 1 
are shown, except those for the MC, which are shown for 
completeness, but without their error bars. In the limit 
of large L, we find for WS bond, MC, and BL data an 
agreement with the range of values for J7 obtained from 
the fit quality method, which corroborates our numerical 
results. For the WS site model in 2D ficst is consistent 
with w = 0.6, while in 3D it agrees with the other mod- 
els. We cross checked also by applying other methods 
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FIG. 6. (color online) Estimated leading correction Hcst as 
defined in Eqs. ([7| and ||8| from the mass data of the bridge 
line, the watershed on bonds (sites), and the main crack in 
2D. The value of a is fixed to 1.2168, the fractal dimension 
estimated by the fit quality method. For better visibility, the 
data of each model is shown with connecting lines and data 
points with Af2ost > 1 have been removed. The values for 
the main crack (MC) are shown for comparison, but without 
their error bars. The horizontal lines give the value (solid) 
and error bar (dashed) for Q, as estimated by the fit quality 
method, as well as the value for u (dotted). 



like, e.g., the one used in Refs. 



and found results 



consistent with the ones presented here. 



IV. RELATION BETWEEN THE MODELS 
A. Bridges, Cracks and Great Walls 

The numerical agreement between the bridges in 
ranked percolation, the optimal path cracks in the strong 
disorder limit, and the watersheds in the 'great wall 
model' supports the claim, made in Sec. |ll| that these 
models are completely equivalent. More precisely, they 
correspond to different strategies for finding the same ob- 
ject (the watershed, the bridge line, and the optimal path 
cracks, respectively). Since these strategies also use the 
random number generators in different ways, they lead to 
different statistical errors, but they give identical scaling 
laws and identical corrections to scaling. 

The random occupation procedure in ranked percola- 
tion !28j can be interpreted as rank sites by increasing 
order in the energy and iteratively occupy them accord- 
ing to their position in the rank. At every step, each 
occupied site has a lower energy than any unoccupied 
one. In strong disorder, the energy of any path is domi- 
nated by the one of the site with the largest energy and, 
therefore, a path of occupied sites, has always lower en- 
ergy than any path containing unoccupied ones. Occupy- 
ing the first bridge site would lead to a spanning cluster 
(SC) and for the first time enable paths that connect the 



FIG. 7. (color online) Estimated leading correction fiest as 
defined in Eqs. ([7| and (|8| from the mass data of the bridge 
line, the watershed on bonds (sites), and the main crack in 
3D. The value of a is fixed to 2.487, the fractal dimension 
estimated by the fit quality method. For better visibility, the 
data of each model is shown with connecting lines and data 
points with Aficst > 1 have been removed. The values for 
the main crack (MC) are shown for comparison, but without 
their error bars. The horizontal lines give the value (solid) 
and error bar (dashed) for Q, as estimated by the fit quality 
method. 



two opposite borders. The bridge site, as being the last 
occupied one in such a path, has the largest energy of 
all sites in it and characterizes the energy of the path. 
The optimal path is one of those paths, as their energy 
is lower than any other connecting path passing through 
unoccupied sites. This means that the first optimal path 
is cracked at the bridge site. Proceeding with the oc- 
cupation of sites following the rank, the next time con- 
necting paths are obtained is when the next bridge site 
is occupied. Again, the energy of the new optimal path 
is dominated by the energy of the current bridge site. 
As before, the crack appears at the bridge site. In this 
picture, the optimal paths always crack at bridge sites, 
until the system is completely disconnected. We, there- 
fore, conjecture that the bridge line and the optimal path 
crack are identical. 



B. Interrelations between the Three Watershed 
Models 

As also seen from the different corrections to scaling, 
the relationships between the three watershed models are 
less trivial and, indeed, quite subtle. 



1. Bond and Great Wall Models 

Both in the bond model and in the great wall model, 
watersheds are topologically strictly one-dimensional 
closed chains. Removing even a single bond (site) from 



them would cut them open, and removing two non- 
adjacent bonds (sites) would cut them into two disjoint 
pieces. Furthermore, one can easily see that any bond in 
the bond watershed must be dual to a bond adjacent to 
a site in the great wall, and that any such site can have 
at most three adjacent bonds corresponding to bonds in 
the bond watershed. This gives immediately 



-Mbond < SMgreatwall, 



(9) 



^/, bond 



< 



and therefore also the rigorous inequality d 

^/, grcatwall- 

We have no similar argument for the opposite inequal- 
ity, but our numerics suggest of course strongly that both 
fractal dimensions are the same. 
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2. The Site Watershed Model 

Although one might have anticipated that the great 
wall model is more similar to the site model than to the 
bond model, the opposite is true. Indeed, the site model 
shows a strong anomaly that makes its finite size cor- 
rections very different, although it seems that it still has 
the same fractal dimension. This anomaly is clearly seen 
in Fig. [8] where we compare the cumulative mass distri- 
butions obtained for BL, WS bond, and WS site of 2D 
systems with size L = 128. While these distributions fall 
off rapidly (roughly Gaussian) for the BL and WS bond 
models, we see a very pronounced tail in case of the WS 
site model. Similar observations have been made, e.g., in 
Ref. [45j. This tail still falls off fast enough to have no 
effect on the fractal dimension, but it definitely calls for 
an explanation. 

Indeed, the watershed in the site model is not strictly 
one-dimensional in the topological sense, but can con- 
tain arbitrarily "thick" regions where it is effectively two- 
dimensional. These regions correspond to lakes with a 
single outlet site, from which the water can run off to- 
wards both sinks. Their existence can also be deduced 
from the hooding algorithm used to construct the site 
model watershed: As explained in Sec. Inj any site 'up- 
stream' of a watershed site has to be also on the water- 
shed. An example of a very small system showing this 
phenomenon is given in Fig. [91 As exemplified in this 
figure, it follows from the algorithm that the great wall 
is always a subset of the site model watershed. Thus one 
has the strict inequalities 



Msitn > M, 



grcatwall ; 



(10) 



and d/, site > c'/, grcatwall ■ Again we cannot prove rigor- 
ously the opposite inequality for the fractal dimensions, 
but again the numerical evidence for equality is over- 
whelming. 

The origin of the power-law tail lies deep in the defi- 
nition of the watershed on sites, namely in the fact that 
entire branches in the diverting runoff scheme can be 
part of the watershed. We will explain this here for the 



FIG. 8. (color online) Cumulative distribution P(ms > m) of 
the masses obtained with the WS site, WS bond, and the BL 
model for system size L = 128 in 2D. The tail of the WS site 
case follows a power law with exponent —1.8. 



representative system depicted in Fig. [9] First, we start 
with the BL, occupying the sites in increasing order of 
the heights, so 1, 2, 3, . . . , 9. The first percolating cluster 
we would obtain when 6 is occupied, which is therefore a 
bridge site and the same applies to 7 and 8, while 9 just 
belongs to the bottom part. For WS site, we find that 
starting from 6 three branches develop, one to 1 and the 
bottom sink, another to 2, passing to 4 and reaching the 
bottom sink (passing 1 and directly from 4) , and a third 
to 3 and the top sink. Hence, from 6 both sinks can be 
reached and it is therefore part of the watershed. The 
same is true for 8, two going to the top (3 and 5) and 
one to the bottom (2). If we now start our runoff scheme 
from 9, we see that initially it diverts into four branches, 
where three are part of the basin of the bottom sink. 
But the branch going upwards, is split at 7 into three 
sub-branches (to 2, 5, and 6), the branch from 2 again 
reaches the bottom sink, but the one growing from 5 is 
part of the top basin. Hence, 7 is part of both (or neither) 
basin, so it is part of the watershed and, by definition, 
also its parent 9 has to be considered part of the water- 
shed. Similarly, this can be deduced from the sub-branch 
to 6. The watershed of this system, therefore, consists of 
the BL and an overhang of one additional site. In gen- 
eral, such overhangs can be larger than one site but all 
bridge sites are always part of the watershed. 

As we have conjectured, the main crack and the bridge 
line are identical, such that discussing the relation of WS 
site and BL is equivalent to discuss the relation of WS site 
and MC. Considering the elevations of a landscape to cor- 
respond to energies, e.g. potential energy, its watershed 
and its optimal path crack can be compared. We have 
defined that a site belongs to the watershed, when the 
invasion percolation clusters grown from two lower near- 
est neighbor sites do reach the opposite borders (catch- 
ments). As both clusters, by definition, do not cross the 
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FIG. 9. Representative system with L — 3, where each square 
cell represents a site of the lattice. The numbers in the lower 
left corner of each cell give the heights. Letters "B" and "W" 
indicate that a site is part of the bridge line (i.e. the great 
wall) and of the site watershed, respectively. Notice that the 
center site in the bottom row is part of the site watershed (as 
it is upstream of the central site), but is not part of the great 
wall, because the wall built at the center site prevents water 
to flow there. Arrows indicate the flow of water from sites 
belonging to the two basins. 



watershed, the watershed site separating the two, has 
a larger height (energy) than any site in both clusters. 
Therefore, for each watershed site (also the overhangs) 
there are always paths which consist solely of sites with 
lower energy than the watershed site, connecting it to 
either border. In the strong disorder limit, the energy 
of these paths is dominated by the largest local energy, 
i.e. the energy of the watershed site. Hence, every path 
connecting the two borders has at least the energy of the 
site where it crosses the watershed and the optimal path, 
the one of lowest energy, crosses the watershed at its low- 
est site. From the same arguments it follows that it is 
then the watershed site which is removed by the OPC 
model. After this, the next optimal path will cross the 
watershed at the next site in increasing order in the en- 
ergy and cracks at the watershed as well. Until the final 
disruption of the system, in strong disorder, every crack 
appears at the watershed site. Hence, the MC is also a 
subset of the watershed. 

From these findings it follows that the power-law tail 
in the mass distribution for WS site arises due to the 
existence overhangs. We know from previous studies, 
that the finite-size cut-off of distributions which follow 
a power law can heavily affect the scaling behavior of 
the moments of this distribution |13j . The average mass 
M is the first moment of the mass distribution P{ms) 
(the derivative of the distribution shown in Fig. ^ and, 
therefore, its scaling behavior is affected by the cut-off 
L^ of its power-law tail. As we based our analysis of the 
corrections to scaling on M, also Cl might be affected. 
We observe the upper cut-off of the tail to scale with 
L? and the lower cut-off to scale with L'*-'. Therefore, 
the functional form of the tail of the cumulative distribu- 
tion is given by P{ms > m) oc m~^'^L^'^'^f . To quantify 
the contributions of the overhangs to C^, we derive here, 
similar as it was done in Ref . jl3] , the contribution of the 
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(11) 



const) , 

what leads to a contribution of order L^'^-^ to the correc- 
tions to scaling of WS site. Although it is only a rough 
estimate, the similarity of this contribution to the value 
we found for the leading correction (w ~ 0.6) is striking. 
The other models have no such overhangs and therefore 
the corresponding amplitude is very small. Together with 
the fact that for these other models the amplitudes of the 
u! correction are small, this suggests that this term in Cl 
of WS site only arises due to the overhangs. Apart from 
this we find the corrections to scaling of all models to 
be in agreement with each other. Furthermore, in 3D no 
such power-law tail is observed for the watershed on sites 
and all models hence have similar distribution of masses. 



V. CONCLUSION 

We obtained from a correction-to-scaling analysis, with 
high precision, an estimate for the fractal dimension of 
the watershed on bonds (WS bond), the watershed on 
sites (WS site), the bridge line (BL), and the main crack 
(MC). We found these fractal dimensions to be, within 
the error bars, in agreement with each other. All models 
have within error bars the same leading correction-to- 
scaling exponent in 2D (second leading exponent for WS 
site) and in 3D. These results are also corroborated by 
the analysis of the local logarithmic slopes in the limit 
of large system sizes. We estimate for all models df = 
1.2168 ± 0.0005 and f] = 0.9 ± 0.1 in two dimensions and 
df = 2.487±0.003 and Q = l.OiO.l in three dimensions. 
The equivalence between the models is also supported 
by either heuristic or exact arguments. Furthermore, we 
give an explanation for the origin of the leading correction 
for WS site in 2D. The estimated values agree with the 
fractal dimensions obtained in previous studies for the 
watershed [TT1I13) . the optimal path crack [14Ul6j . and 
the bridge line p5H5D] . as well as with the ones found 
for the perimeter of the infinite cluster in discontinuous 
percolation (1.23 ±0.03 [34 and 2.5 ±0.2 [35 ). It would 
be interesting to know if this perimeter also obeys the 
same corrections to scaling as we have found. 
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